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Abstract —Function optimization and finding simulta- 
neons solutions of a system of nonlinear eqnations (SNE) 
are two closely related and Important optimization prob¬ 
lems. However, unlike in the case of function optimization 
in which one is required to find the global minimum 
and sometimes local minima, a database of challenging 
SNEs where one is required to find stationary points 
(extrama and saddle points) is not readily available. 
In this article, we initiate building such a database of 
important SNE (which also includes related function 
optimization problems), arising from Science, Engineering 
and Economics. After providing a short review of the 
most commonly used mathematical and computational 
approaches to find solntlons of such systems, we provide 
a preliminary list of challenging problems by writing 
the Mathematical formulation down, briefly explaning the 
origin and importance of the problem and giving a short 
acconnt on the currently known resnits, for each of the 
problems. We anticipate that this database will not only 
help benchmarking novel numerical methods for solving 
SNEs and function optimization problems but also will 
help advancing the corresponding research areas. 

L Introduction 

Development of methods for local and global op¬ 
timization, which include finding the global minimum, 
local minima and saddle points, of nonlinear multivari¬ 
ate objective functions, say F{x), has always been one 
of the most active areas of research in Mathematics and 
Computer Sciences, due to their applications in many 
areas of Science, Engineering, Economics, etc. Here, 
F(x) is usually a real-valued function from to R. 
The A^-dimensional space is made of the degrees of 
freedom of the physical system. 


The most general form of the optimization problems 
is to find the stationary points (SPs) of F(x), defined as 
the simultaneous solutions of the system of equations 
/i(x) = dF{x)/dxi = 0, for all i = The 

SPs at which exactly i eigenvalues of the Hessian are 
negative definite, and the remaining N — i eigenvalues 
are positive definite, are called saddles of index i, with 
* = 0 SPs also known as the minima of E'(x). The 
SP at which F{x.) attains its lowest value is known as 
the global minimum, provided F(x) is bounded from 
below. The SPs with at least one 0 eigenvalue are called 
singular SPs. 

As is common to most nonlinear problems, an 
analytic calculation of the SPs is extremely difficult, 
and in most cases impossible. Hence, one has to rely 
on computational methods. Leaving the certification of 
numerical solutions aside Q-ii, the Newton-Raphson 
(NR) approach has been a popular method to find SPs of 
nonlinear functions. There, one refines an initial guess 
through successive iterations and hopes to converge 
to a solution. However, the NR method may often 
converge slowly, or even worse but not uncommon, may 
diverge depending on the initial guess, and may behave 
erratically near singular solutions a, a. 

One can also resort to alternative methods such as 
the gradient-square minimization method: instead of 
solving /i(x) = 0 directly, here one minimizes the 
sum of squares W = using traditional 

numerical methods, such as conjugate gradient a, 
0. When imposing the further constraint W = 0, 
the corresponding minima of W are then the desired 
solutions of /i(x) = 0, for all * = 1,..., iV. However, 
the number of minima with W > 0, which are not 
the solutions of /i(x), generally outweighs the desired 
minima at which W = 0. Moroever, these non-solutions 


may also be singular making the minimization problem 
ill-conditioned ||9l, ifTOl . Hence, this approach turns out 
to be very inefficient in practice a, CD. Instead, a 
biased gradient squared descent framework ifT^ may 
provide a more useful alternative to the gradient squared 
minimization method. 

There are other systematic approaches such as the 
one based on the Broyden-Fletcher-Goldfarb-Shanno 
(BFGS) algorithm CS, ifTD . or eigenvector-following 
a, CD (implemented in the OP TIM package, which 
also includes many other geometry optimization tech¬ 
niques, such as a modified version of the limited- 
memory BFGS algorithm ca, d, single- and double- 
ended ca-na searches). 

For SNEs having polynomial-like nonlinearity, sym¬ 
bolic algebraic geometry methods based on the Grobner 
basis technique can guarantee to find all the complex so¬ 
lutions ll^ . However, the algorithms are known to have 
exponential space complexity, which may make solving 
even a moderately sized SNE prohibitively difficult. 
Another rigorous approach which can guarantee to find 
all the solutions of a system of nonlinear equations is 
an interval based method f22\ . However, this approach 
has only proved successful for a very small systems 
and SPs so far because it is based on bisections of 
the ranges and the computation blows up by increasing 
the number of variables. In the homotopy continuation 
approach ED, on the other hand, one starts with a 
new SNE that is qualitatively similar to the SNE to be 
solved, and whose solutions are known or can be easily 
obtained. Then, each solution of the new system is 
homotopically continued to eventually obtain a solution 
of the system to be solved. Eor a general nonlinear 
system, this method does not usually guarantee to find 
all solutions Q. 

Recently, a specialized homotopy continuation 
method based on algebraic geometry, namely the numer¬ 
ical polynomial homotopy continuation (NPHC) method 
Il2^ - ll28l . has captured the attention due to its ability of 
finding all the isolated SPs of F{x.) having polynomial¬ 
like nonlinearity. Several good implementations of the 
NPHC method are now available ED . E9ll - BT]| and the 
method is applied to many different areas in Science and 
Engineering in recent years EDl - ll44| . In this method, 
after coming up with an upper bound on the number of 
isolated complex solutions of the given SNE, the system 
is continuously deformed from a different system whose 
solution count agrees with the upper bound to finally ob¬ 
tain all the complex solutions of the original system. In 
the end, only real solutions are retained being physically 
relevant. However, in many real life applications, the 
number of complex solutions may be extremely large 
making the task of computing all of the real solutions 


'There have been attempts to construct specialized homotopies 
which guarantees to find all isolated solutions for general (i.e., non¬ 
polynomial) SNE (see, e.g.. l241 . 1251 ). though these claims have not 
yet been rigorously proven. 


a prohibitively difficult task. 

An inversion-relaxation method El was also re¬ 
cently introduced in which, to obtain stationary points 
of a given F{x.) which is bounded from above and 
below, one relaxes from a solution of index N to find 
all the saddles of index {N — 1) that are connected 
to the maximum. Relaxing from each of these saddles 
of {N — 1) index, one then obtains all the saddles of 
{N — 2) index that are connected to the corresponding 
saddle of (A^ — 1) index, and so on. One can then obtain 
many saddles of all the possible indices starting from 
one maximum. If all the maxima of F(x) are known, 
then all the stationary points may be found using this 
method. Eurther investigations in addition to a rigorous 
proof of the previous statement is still needed. 

In recent years, soft computing methods (especially 
population based methods) prove to be efficient in find¬ 
ing multiple solutions with relatively less information 
about the system (i.e. without using derivatives, etc.). 
But there are still challenges for these methods, par¬ 
ticularly in the following aspects: dealing with systems 
which have a large number of equations; finding a large 
number of distinct solutions; and scaling an efficient 
method for a simpler system to larger systems. The work 
in ll46l is among the first approaches which transform 
a system of equations into a multiobjective optimiza¬ 
tion problem. Although various ways to transform a 
system of equations into an optimization problem have 
been proposed El-lED, none of these methods can 
detect all the solutions, and are biased for problems 
in specific areas. The main ways in which population 
based heuristics transform a system of equations into an 
optimization problem are single-objective optimization 
based methods, constrained optimization based meth¬ 
ods, and multiobjective optimization based methods. 

However, as with many other Mathematical areas, 
optimization method development is mostly discon¬ 
nected to the real life applications. Most of the con¬ 
ventional test systems 1581 currently used by the opti¬ 
mization community, though while serving the purpose 
of benchmarking the novel methods may not appeal 
the Scientific community. Another reason for the dis¬ 
jointedness is the different terminology, related to the 
corresponding scientific application, being used among 
the scientific communities for the same underlying 
Mathematics. 

Our goal in this paper is to bridge the gap by col¬ 
lecting some of the challenging systems of equations as 
optimization problems arising from real life applications 
in the terminology which is accessible to the optimiza¬ 
tion and computational mathematics community. Hence, 
our list of challenging systems may not only provide 
motivation and benchmarking for novel algorithms, but 
solving them for the yet unsolved cases will provide 
advancement in the respective areas the systems arise 
from. The models we provide are described generic, 
given by their general representation. They can be used 



as benchmarks of variable size, some being simpler 
(lower dimension) than the others. All the systems 
have multiple solutions, in some cases this increases 
exponentially with the dimensionality of the system. 
Finding all solutions will be one of the main challenges 
of the population based algorithms. As such a set of 
benchmarks is currently missing in the evolutionary 
computation literature, we believe that our work herein 
will be of help for all the researchers working on 
computational approaches (especially in evolutionary 
computation field) for solving complex systems of equa¬ 
tions. 

In the remainder of this paper, we provide a brief in¬ 
troduction, Mathematical formulation, and a list of open 
Mathematical problems for each of the optimization 
models We do not however intend to provide a complete 
list of references nor a complete list of available results 
for each of the problems. Rather we refer the reader to 
the corresponding databases, if available. 


II. Challenging Problems 

In this Section, we list out a few of the challenging 
problems arising in various research areas. 


A. Chemical and Physical Clusters 

In theoretical chemistry, physics and many other 
areas in Science and Engineering, exploring the hyper¬ 
surface defined by a multivariate function, V (x), called 
the potential energy function, plays a very important 
role in understanding and describing the physics and 
chemistry of the phenomenon. The hypersurface is 
called the potential energy landscape (PEL). 

In fact, a variety of methods based on the SPs have 
recently attraced a lot of attention of both chemists 
and physicists, due to their applications to many-body 
systems as diverse as metallic clusters, biomolecules, 
structural glass formers, and coarse-grained models of 
soft matter, etc. ll59l , IbO). Einding SPs of V (x) provides 
the foundations for global optimization llbTI - ll^ , ther¬ 
modynamic sampling to overcome broken ergodicity 
JMI-ISI, as well as rare event dynamics Esi-ca 
within the general framework of PEL theory ||59|. 
Below, we list out a few important potentials coming 
from chemistry and physics applications as optimization 
problems. 

1) The Nearest-neighbour Model: The two- 
dimensional nearest-neighbor (j)^ model has been widely 
studied because (1) it is one of the simplest models with 
a continuous configuration space, (2) it exhibits a phase 
transition in the same universality class as the two- 
dimensional Ising model. Eor an N € and J,X,^ G 
R the model, in variables x = (xii,xi 2 , ■ ■ ■ ,xnn), 


is V (x) given by 


\ 

[xij Xf^l ) 

( 1 ) 

where A C is the standard square lattice with N'^ 
lattice-sites and C A is the four nearest-neighbor 

sites of The stationary equations are given by 
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for each pair of i,j = 1,... ,N. The traditional bound¬ 
ary conditions is the periodic one, A = 3/5 and /i^ = 2. 
Only real solutions are physical for this model. 


The model has played a crucial role in studying link 
between the topology of the potential energy landscape 
with the phase transition ll^ . II^ . ITT^ . Moreover, the 
model shows an interesting behavour while varying the 
parameter J from 0 to 1, i.e., one can go from the case 
when all the solutions are real to only an extremely 
small fraction are real. A variety of computational 
methods have been used to explore the PEL of this 
model. 


The NPHC method has found all the SPs for N = 
3,4 in a previous study m, d- However, this model 
poses a particularly tough challenge to the method since 
the total number of solutions in , counting multi¬ 
plicity, is always equal to its total degree (the Bezout 
bound) 3^ , which grows quickly as N increases. Eor 
example, for N = 6 and 7, the total degree exceeds 
10^^ and 10^^, respectively. Hence, finding all complex 
solutions is clearly unfeasible with current technology 
for large values of N. More recently, the Newton 
homotopy method has been employed to find many real 
solutions of this model for larger N, though without 
a guarantee of finding all the real solutions ||6l. Thus, 
this model still needs to be studied with more refined 
methods. On the other hand, due to the simplicity of 
the model, it can be used as a benchmark model. 

2) The Thomson Problem: The Thomson problem 
ca is to find the minimum energy of the system 
made of N electrons restricted to move on the surface 
of a sphere of unit radius. The model was originally 
proposed by Sir JJ Thomson as a natural consequence 
of his atomic model called the plum-pudding model. 
Though the model turned out not to be a correct model 
for the atoms due to experimental evidence, it has turned 
out to be very interesting model in chemistry, physics 
and mathematics. The global minima found by numer¬ 
ical methods have been observed to be geometrically 
irregular: though the global minima for = 4, 6 and 
12 are the expected platonic solids, surprisingly, those 
for N = 8 and 20 are not platonic solids. The Thomson 
problem and its generalizations have also been used to 
model clusters of proteins on a shell, colloid particles. 





fullerene patterns of carbon atoms, etc. In addition, the 
model has got the status of a standard benchmark system 
for any new optimization routines. Finding the exact 
global minimum of the model has also been an active 
area of research, though the exact results are known 
only for a handful of TVs (TV = 2 - 6 and 12')i76l. liTTl. 

Because the electrons interact with each other with 
the Coulomb potential, the potential energy function of 
this system is given by 


one of the authors has worked out a Newton homotopy 
method which can not only efficiently find multiple 
stationary points of this model, but also can find singular 
solutions of the model. The latter task has been proved 
to be prohibitively difficult using traditional methods 
in which one has to invert the hessian matrix at such 
singular solutions. 

4) Morse Clusters: The Morse potential ||8^ is 
given as 


^Th(^)= E (3) 

l<i<j<N 

where = y'{xi - XjY + (t/i - yjY + (z* - ZjY, 
with constraints xf +yf + zf = 1 for alH = 1 ,... ^N. 
To remove the rotational symmetry of the system, we 
fix xi = = 2/2 = 0 and zi = 1. 

The model has gained special attention recently as 
it appeared as the 7th problem in Steven Smale’s list 
of eighteen unsolved problems for the 21st century 
GH). This problem was indeed motivated by finding 
a good starting polynomial for a homotopy algorithm 
for realizing the Fundamental Theorem of Algebra. For 
the sepcial case of finding the global minimum of the 
model, in Ref. Il79l it was shown that good initial points 
to find the global minimum of the Thomson model are 
in fact complex roots (projected to 2-sphere) of random 
univariate polynomials. 


3) Lennard-Jones Clusters: One of the most popular 
model for atomic interactions, in theoretical chemistry 
is the Lennard-Jones potential |[80l, which is defined as 



where e is the pair well depth, and 2 ^/®ct is the 
equilibrium pair separation. We take e = a = 
1. Here, same as in the Thomson problem, ry = 
yj{xi — Xj)"^ + (t/i — UjY + {zi — ZjY is the distance 
between atoms i and j. To remove the global degrees 
of freedom of the model coming via the rotational 
and translational invariance of the system, we can fix 
xi = yi = zi = y 2 = Z 2 = Z 3 = 0. Thus, in total there 
are 3iV — 6 variables in l^j yielding 3iV — 6 equations 
Vl^j = 0. The model has gained popularity among 
the theoretical and computational chemists because it 
is simple enough to perform molecular dynamics sim¬ 
ulations and also because it is fairly accurate approx¬ 
imation of the actual atomic interactions validated by 
experiemental observations. The minima of this model 
represent individual molecular configurations. There is 
a huge literature available addressing the global and 
local optimization issues of the PEL of this model. An 
extensive search for minima and saddle points has been 
carried out in © for N up to 13, and a search for 
minima and saddles of index one (transition states) for 
N = 14 was presented in lISTl . Recently, in Ref. Il82l . 


F(r) = -2)), (5) 

i<j 

where ry is again the distance between atoms i and j, e 
is the dimer well depth and equilibrium bond length. 
Since these two parameters do not affect the geometry 
of the PEL, we can conveniently set them to be unity. 
We can again fix 6 of the coordinates to remove the 
rotational and translational symmetries of the system 
to have eventually a system of 2>N — 6 equations and 
‘3N — 6 variables. 

The most important parameter in this model is p 
which is dimensionless, p determines the range of the 
inter-particle forces: low values of p correspond to long 
range interactions because increases the range of the 
attractive part of the potential and softens the repulsive 
wall, thus widening the potential well. Similarly, 
large values correspond to short range interactions. 
The ability to continuously varying the interaction 
range of the particles has made this model widely 
popular to study a range of chemical phenomena from 
intermolecular potential of Cqq to alkali metals. In fact, 
at p = 6, the above potential has the same curvature at 
the bottom of the well as the Lennard-Jones potential. 
A database of the known minima and transition 
states of this model for various values of p are 
available at the Cambridge Energy Landscape database 
http://www-wales.ch.cam.ac. uk/ CCD.html 

5) The XY Model: The XY model is one of the 
simplest potentials with continuous degrees of freedom 
(unlike the Ising model in which the configuration space 
is discrete) in theoretical physics and chemistry, though 
its potential energy landscape is very rich and inter¬ 
esting, and has been helpful in understanding general 
features of potential energy landscapes. 

The function to be extremized for the XY model 
on a d-dimensional cubic lattices A of side length L is 

d? = ^ E E cos(6lfc -9i)]. (6) 

fceA i&M(k) 

Here, the total number of lattice sites is N = L‘^, 
and for each lattice site k € A the angular variable 
6 k € (—TT, tt]. denotes the set of nearest-neighbors 

of lattice site k. Moreover, the parameters Jk^i are 
the random disorders which are i. i. d. picked from 
some random distribution. One has to pick a boundary 
condition here due to the nearest-neighbour terms, the 






usual choices being periodic and anti-periodic boundary 
conditions. For the periodic boundary conditions, to 
remove the global 0(2) symmetry, exactly one of the 
angles is fixed to zero. 

The function (ISll also appears in many other areas 
such as in statistical physics ll84ll . complex systems ll85l . 
lattice field theories ll^ . [861, IISTl . etc. The model 
is used to model low-temperature superconductivity, 
superfluid helium, hexatic liquid crystals, and other 
phenomena. 

For the one-dimensional case, all SPs of this model 
are analytically found in 1871 - 1891 . Using the SPs 
of the one-dimensional model, a class of SPs for the 
two- and three-dimensional cases can be constructed 
mni (see also (Ml, (HI)- For the two-dimensional 
case, for small number of lattice-sites, all the SPs were 
found using the NPHC method in l32l. l38l. l87l. l92l. 
Using other traditional numerical methods (Ml, dQi, 
it was then shown that the number of isolated SPs, 
as well as the number of minima, of the XY model 
in two- and three-dimensions increases exponentially 
as N increases. Similar results were obtained from 
the Kuramoto model point of view in l93l . It was 
also shown that even after removing the global 0(2) 
symmetry, the model possess many continuous SPs. 
Several attempts for finding the global minimum of 
this model have also been made, e.g., using Simulated 
Annealing in (941 . However, the model has proven to 
be a very challenging optimization problem and even 
for moderate N many features of the PEL is yet to be 
explored. 

B. Polynomial systems in Economics 

The computation of equilibria in economics leads 
to systems of polynomial equations. Also known as 
pseudo-games, social equilibrium problem, equilibrium 
programming, coupled constrained equilibrium prob¬ 
lem, abstract economy (951 . Using the notations and 
definitions from (9^ . the n-person game is defined by 
n players, labelled 1, 2, ..., n. Each player i has di pure 
strategies labeled 1, 2, . . . , d,;. Each player has an ob¬ 
jective function that depends both on his own strategies 
and the strategies of the other players. This function is 
called utility function or payoff function. The game is 
defined by n payoff matrices ...,X^^\ one 

for each player. Each matrix X^'^'i is an n-dimensional 
matrix of format dl x d2 x ... x dn whose elements are 
rational numbers. The element represents the 

payoff for player i if the other players 1, 2,..., n select 
strategies respectively. Each player selects 

a mixed strategy given by 

where is the probability of player i to select strategy 
j. The vector p*^®^ is a probability distribution of player 
i on his set of pure strategies. 

The payoff tt^ for player i is given by his matrix 

A:(®): 


dl d 2 

E ■ 

ii=i j2=i 


^ ^ Pjl Pj2 "'Pjn 

jn = l 


Jn) 


Since the variables of the problem are probabilities, 
we have that p^-®^ > 0,Vi, j, and p^®^ + P 2 + = 

1,Vj which means that p = (p*-®^) is a point in the 
product of simplices A = Adi — lx Ad 2 — 1 x ... x 
Ad„ - 1. 


A point p S A is a Nash equilibrium if none of the 
players can increase his payoff by changing his strategy 
while the other n — 1 players keep their strategies 
unchanged. This can be expressed as a system of poly¬ 
nomial constraints in the variable vectors p S A and 
TT = (tti, ..., 7r„) £ S'®, with the following multilinear 
polynomial for each p^®^: 


Pk^i 


dl 

E- 


di-i di+i 

• E E ■ 

ii-i=iii+i=i 


dn 

Jlj2-jn 

3r. = l (7) 


X„(1) X 77^®+^^ n^®®^ 

^Pjl ■■■Pji-l ^ Pji + 1 --Pjn 


which, together with the constraints, represents a system 
of n-|-di-|-...-|-d„ equations in n-|-di-|-...-|-d„ variable. 
Each polynomial is the product of a linear polynomial 
and a multilinear polynomial of degree n — 1. 

A solution (p, tt) £ A x 5R®® represents a Nash 
equilibrium for the game if and only if (p, tt) is a zero 
of the polynomials in (Tji and each expression in the 
parenthesis is nonnegative. 


C. Edge matching puzzles 

Edge-matching problems are popular puzzles in 
which, given a set of pieces and a grid, the goal is 
to place the pieces on the grid such that the edges 
of the connected pieces match. Edge-matching puzzles 
are challenging because there is no global image as 
guidance and there is no guarantee that two pieces 
fitting together are in the right positions. Edge-matching 
problems are proved to be NP-complete (971 . A direct 
implication and practical application of the edge match¬ 
ing puzzles is in image reconstruction. 

The work in (981 formulates the edge matching 
problem as a systems of polynomial equations derived 
from the pieces of the puzzle. Solutions of the system 
represent solutions of the puzzle. The authors consider 
a particular instance of the problem in which a set of 
pieces of known shapes and edge colours, is given. The 
puzzle is bounded by a frame and each edge must match 
either an edge of another piece or the frame. The puzzle 
has N pieces and the puzzles as considered as a 2- 
dimensional polygons. Each piece i of the puzzle is 
represented by its location ti £ given by its center 
and its set Ei of edges. Each edge j £ Ei is described 
by the relative location of its center bij with respect 
to the piece center, its color Cij and its orientation 
or inclination given by the angle 6 ij . The absolute 


location of the jth edge of the ith piece is given by the 
sum ti + bij. 

The pieces corresponding to the puzzle frame have i 
= 0 and the properties of its edge elements 6oj, Cqj and 
Oqj. The only operation that can be applied to puzzle’s 
pieces is translation. The goal of the game is to find 
a translation such as all edge elements pair 

with matching edge elements in their spatial location, 
colour and orientation. 

If (fi,...,fAr) is a solution of the puzzle then it is 
also a solution for the system of equations given by: 

for every (c, 9) and every real valued function / : 

5R, where Sij{c,9) is a signed indicator function w.r.t. 
(c, 9) and is given by: 

{ 1 , ^ 

-1, Cij = c, = 6> + TT (9) 
0, otherwise 

Different choices of /-function are possible. There 
trivial one reduces the system to: 

'^Sij{c,9){ti + hij) = Q ( 10 ) 

ij 

The fact is that the converse of this statement does 
not hold and not every solution of the system is a 
solution of the puzzle. In order to have the converse 
valid as well, consider / as an exponential function 
defined as (for a given k € ?R.^) fk{u) = In this 

case the equation ([Sll becomes: 

=0 ( 11 ) 

which by replacing 

^s,,,{c,9)Tl^ = 0 (12) 

This leads to the conclusion that fi, is identified 

with Ti,...,Tiy and the fact that the converse holds. 
This, a solution of the system of equation is a solution 
of the puzzle. 

III. Conclusion 

We presented several models which compose com¬ 
plex systems of equations, together with their descrip¬ 
tion. These models arrive in a variety of scientific 
areas and are of great importance. The equations are 
presented in a general format so that those interested in 
using them as benchmarks can build their own instance 
of the system. The scope of this work is to offer 
the research community (especially the computational 
science community) a common set of difficult equations 


systems benchmarks. The systems presented are difficult 
both in terms of number of equations they contain as 
well as in number of alternative solutions. As such a 
repository is mandatory for testing the performance of 
the newly proposed methods for solving systems of 
equations, we believe that the work here represents a 
starting point and that more system models will be 
considered in the future. 
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